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Abstract 

We present a program that implements the OPP reduction method to extract the 
coefficients of the one-loop scalar integrals from a user defined (sub)-amplitude or 
Feynman Diagram, as well as the rational terms coming from the 4-dimensional part 
of the numerator. The rational pieces coming from the e-dimensional part of the 
numerator are treated as an external input, and can be computed with the help of 
dedicated tree-level like Feynman rules. 

Possible numerical instabilities are dealt with the help of arbitrary precision 
routines, that activate only when needed. 
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1 Introduction 



Developing efficient tools to compute one-loop corrections for multi-particle processes is 
an important task needed to cope with the complexity of LHC and ILC Physics. In 
the last few years a big effort has been devoted by several authors to this problem [1]. 
The used techniques range from analytic methods to purely numeric ones, also including 
semi-numerical approaches. 

In the analytical approaches, computer algebra is used to reduce generic one-loop 
integrals into a minimal set of scalar integrals and remaining pieces (called rational terms) , 
mainly by tensor reduction [2-5]. For multi-particle processes this method becomes quite 
cumbersome because of the large number of generated terms. 

In the numerical or semi-numerical methods a direct computation of the tensor in- 
tegrals is performed [6], capable, in principle, to deal with any configuration of masses. 
However, their applicability remains limited due to the high demand of computational 
resources and the non-existence of an efficient automation. 

In a different approach, called the unitarity cut method [7], the one-loop amplitude 
rather than the individual integrals are evaluated, avoiding the computation of Feynman 
diagrams. In another development, the 4-dimensional unitarity cut method has been used 
for the calculation of QCD amplitudes [8], using twistor-based approaches [9]. Moreover, 
a generalization of the the unitarity cut method in d dimensions, has been pursued re- 
cently [10]. Nevertheless, in practice, only the part of the amplitude proportional to the 
loop scalar functions can be obtained straightforwardly. The remaining rational part, 
should then be reconstructed either by using a direct computation based on Feynman 
diagrams [11-13] or by using a bootstrap approach [14]. Furthermore the complexity of 
the calculation increases away from massless theories. 

In two recent papers [15,16], we proposed a reduction technique (OPP) for arbitrary 
one-loop sub-amplitudes at the integrand level [17] by exploiting numerically the set of 
kinematical equations for the integration momentum, that extend the quadruple, triple 
and double cuts used in the unitarity-cut method. The method requires a minimal infor- 
mation about the form of the one-loop (sub-) amplitude and therefore it is well suited for 
a numerical implementation. The method works for any set of internal and/or external 
masses, so that one is able to study the full electroweak model, without being limited to 
massless theories. In [18] the OPP method has been used, in the framework of the unitarity 
cut technique, to explicitly compute the subtraction terms needed not to double count 
the contribution of the various scalar integrals. 

In this paper, we describe a F0RTRAN90 implementation of the OPP algorithm. In 
section 2, we recall the basics of the method and present our solution to compute Rational 
Terms and to deal with numerical inaccuracies. In section 3 we outline the conventions 
used in the program. In section 4 we describe the F0RTRAN90 code that implements 
the method and, in the last section, we discuss our conclusions. Finally, two appendices 
integrate the content of the paper. 
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2 Theory and general features 



2.1 The OPP method 

The starting point of the OPP reduction method is the general expression for the integrand 
of a generic m-point one-loop (sub-) amplitude [15] 

A ® = f, ' A = (?+ft) 2 p ? . (1) 

In the previous equation, we use a bar to denote objects living in n = 4 + e dimensions, 
and q 2 = q 2 + q 2 , where q 2 is e-dimensional and (q ■ q) = 0. N(q) is the 4-dimensional 
part of the numerator function of the amplitude. If needed, the e-dimensional part of 
the numerator should be treated separately, as explained in [19]. N(q) depends on the 
4-dimensional denominators = (q + Pi) 2 — m 2 as follows 

m—1 m—1 

N (q) = E [diiohizia) + d(q; ^0*1*2*3)] II A 

m—1 m—1 

+ Yl [c(io*i«2) + c(g;ioii»2)] II A 

?0<ii<i2 i^io,h,i2 
m—1 m—1 

+ E [&(»oii) + 6(g;ioii)] n A 

«0<U i^io,ii 
m—1 m—1 

+ E [a(«o) + a(5;io)] II A 

m—1 

+ P(?)IIA. (2) 

i 

Inserted back in Eq. (1), this expression simply states the multi-pole nature of any m- 
point one-loop amplitude, that, clearly, contains a pole for any propagator in the loop, 
thus one has terms ranging from 1 torn poles. Notice that the term with no poles, namely 
that one proportional to P(q) is polynomial and vanishes upon integration in dimensional 
regularization; therefore does not contribute to the amplitude, as it should be. The 
coefficients of the poles can be further split in two pieces. A piece that still depend on q 
(the terms d, c, 6, a), that vanishes upon integration, and a piece that do not depend on q 
(the terms d, c, b, a). Such a separation is always possible, as shown in [15], and, with this 
choice, the latter set of coefficients is therefore immediately interpretable as the ensemble 
of the coefficients of all possible 4, 3, 2, 1-point one-loop functions contributing to the 
amplitude. 

Once Eq. (2) is established, the task of computing the one-loop amplitude is then 
reduced to the algebraical problem of fitting the coefficients d, c, b, a by evaluating the 
function N(q) a sufficient number of times, at different values of q, and then inverting 
the system. That can be achieved quite efficiently by singling out particular choices of q 
such that, systematically, 4, 3, 2 or 1 among all possible denominators Di vanishes. Then 
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the system of equations is solved iteratively. First one determines all possible 4-point 
functions, then the 3-point functions and so on. For example, calling q$ the 2 (in general 
complex) solutions for which 

D = D 1 = D 2 = D 3 = , (3) 

(there are 2 solutions because of the quadratic nature of the propagators) and since the 
functional form of d(q; 0123) is known, one directly finds the coefficient of the box diagram 
containing the above 4 denominators through the two simple equations 

N(q±) = [d(0123) + d(g±;0123)] JJ D Mo) ■ (4) 

#0,1,2,3 

This algorithm also works in the case of complex denominators, namely with complex 
masses. Notice that the described procedure can be performed at the amplitude level. 
One does not need to repeat the work for all Feynman diagrams, provided their sum is 
known: we just suppose to be able to compute N(q) numerically. 

The modifications one has to apply to the method when working in d = 4+e dimensions 
are described in the next subsection. 

As a further remark notice that, since the terms d, c, b, a still depend on q, also the 
separation among terms in Eq. (2) is somehow arbitrary. Terms containing a different 
numbers of denominators can be shifted from one piece to the other in Eq. (2), by relaxing 
the requirement that the integral over the terms containing q vanishes. This fact provides 
an handle to cure numerical instabilities occurring at exceptional phase-space points. In 
CutTools such a mechanism is implemented for the 2-point part of the amplitude, as 
described in subsection 2.3 . 

2.2 The rational terms 

The described procedure works in 4 dimensions. However, even when starting from a 
perfectly finite tensor integral, the tensor reduction may eventually lead to integrals that 
need to be regularized 1 . Such tensors are finite, but tensor reduction iteratively leads to 
rank m m-point tensors with 1 < m < 5, that are ultraviolet divergent when m < 4. For 
this reason, we introduced, in Eq. (1), the ci-dimensional denominators A, that differs by 
an amount q 2 from their 4-dimensional counterparts 

A = A + f ■ (5) 

The result of this is a mismatch in the cancellation of the ci-dimensional denominators 
of Eq. (1) with the 4-dimensional ones of Eq. (2). The rational part of the amplitude, 
called Ri [20], comes from such a lack of cancellation and is computed automatically in 
CutTools. 

1 We use dimensional regularization as a regulator. 
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A different source of Rational Terms, called R2, can also be generated from the e- 
dimensional part of N(q) (that is missing in Eq. (1)), and should be added on the top of 
CutTools's results. R 2 can be easily computed by using dedicated tree-level like Feynman 
rules, as explained in detail in [20]. The user's conceptual effort required to provide R2 is 
the same needed to supply the input function N(q). We therefore consider the problem 
of computing R 2 completely trivial and solved once for all. 

The Rational Terms R\ are generated by the following extra integrals, introduced 
in [15] 



/ 



d n q 



q z VK £ 



DiDj 2 



2 , 2 (Pi~ PjY 
m i + m j 



+ 0(e) 



J 



d n q - - = -— + 0(e) 

DiD,,D k Di 6 1 ; 



(6) 



The coefficients of the above integrals are computed in CutTools by looking at the implicit 
mass dependence (namely reconstructing the q 2 dependence) in the coefficients d, c, b of 
the one-loop functions, once q 2 is reintroduced through the mass shift 

2 2~2 

~> m i -Q ■ (7) 

One gets 

b(ij;q 2 ) = b(ij) + qW\ij) , 
c{tjk;q 2 ) = c{tjk)+q 2 c^{tjk). (8) 

Furthermore, by using Eq. (7), the first line of Eq. (2) becomes 

m— 1 m— 1 

V (m \q,q 2 )= Yl [^oV 2 « 3 ;9 2 )+4¥i¥3;? 2 )] II A, (9) 
and the following expansion holds 

m 

V^(q,q 2 )=j:q {2j - A) d^- 4 \q), (10) 

where the last coefficient is independent on q 

S 2m - A \q) =S 2m ~ 4) . (11) 

In practice, once the 4-dimensional coefficients have been determined, CutTools redoes 
the fits for different values of q 2 , in order to determine b^ 2 \ij), S- 2 \ijk) and d^ 2m ~^ . Such 
three quantities are the coefficients of the three extra scalar integrals listed in Eq. (6), 
respectively. 
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A different way of computing d^ 2m ~^ is implemented in CutTools when the Logical 
variable inf is set to .true, in subroutine dp_get_coef f icients and subroutine 
mp_get_coef f icients. In this case the code computes 

= lim g^lM!) (12) 

This limit is numerically quite stable and the computation faster. However, the default 
for inf is .false . . 



2.3 Dealing with numerical inaccuracies 

During the fitting procedure to determine the coefficients, numerical inaccuracies may 
occur due to 

1) appearance of Gram determinants in the solutions for which 4, 3, 2 or 1 denominators 
vanish; 

2) vanishing of some of the remaining denominators, when computed at a given solu- 
tion; 

3) instabilities occurring when solving systems of linear equations; 

In principle, each of these three sources of instabilities can be cured by performing a proper 
expansion around the problematic Phase-Space point 2 . An attempt in this direction is 
described in [16]. However, this often results in a huge amount of work that, in addition, 
spoils the generality of the algorithm. Furthermore, one is anyway left with the problem of 
choosing a separation criterion to identify the region where applying the proper expansion 
rather than the general algorithm. 

The solution implemented in CutTools is, instead, of a purely numerical nature and 
relies on a unique feature of the OPP method: the fact that the reduction is performed 
at the integrand level. In detail, the OPP reduction is obtained when, as in Eq. (2), the 
numerator function N(q) is rewritten in terms of denominators. Therefore N(q) computed 
for some arbitrary value of q by using the 1. h. s. of Eq. (2) should always be numerically 
equal to the result obtained by using the expansion in the r. h. s. This is a very stringent 
test that is applied in CutTools for any Phase-Space point 3 . When, in an exceptional 
Phase-Space point, these two numbers differ more than a user defined quantity (limit), 
the coefficients of the loop functions for that particular point are recomputed by using 
multi-precision routines (with up to 2000 digits) contained in CutTools [21]. The only 
price to be payed by the user is writing, beside the normal ones (namely written in double- 
precision), a multi-precision version of the routines computing N(q), that is anyway easily 



2 From now on we will denote such a point as exceptional. 

3 The arbitrary, complex, 4- vector q used for this test is randomly chosen by the code in a point by 
point basis. 
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obtained by just changing the definition of the variables used in the routines, as explained 
in appendix A. The described procedure ensures that the coefficients of the scalar loop 
functions are computed with the precision given by limit. This is usually sufficient; 
however, when strong cancellations are expected among different loop functions, a multi- 
precision version of the one-loop scalar functions should also be used. Then, a complete 
control over any kind of numerical inaccuracy is guaranteed. Finally, one should mention 
that, usually, only very few points are potentially dangerous, namely exceptional, so that a 
limited fraction of additional CPU time is used to cure the numerical instabilities, therefore 
compensating the fact that the multi-precision routines are by far much slower than the 
normal ones. This procedure has been shown to work rather well in practice. 

A final remark is in order. For strictly massless momenta, all Phase-Space points are 
exceptional in the 2-point sector. Differently stated, expressing tensors such as 

or l^fk (13) 

in terms of scalar 2 and 1-point functions necessarily involves the appearance of powers 
of ^ J po )2 ; that is always a problem when (pi — p Q ) 2 = 0. 

For this reason, a different basis [16] is implemented, for the 2-point sector, in CutTools. 
This basis makes use of an arbitrary massless vector v and the code computes the coeffi- 
cients of the following three scalar integrals 

rfra - [(g + Po) - V Y with £ = 0,1,2 and v 2 = . (14) 
DqD\ 

Notice that, when k\ = (pi — p ) 2 = and m Q = mi, 



I " * D Q D l 2 J " H D Q D l ' 



exactly. 



3 Conventions used in the program 

The information to be provided by the user is 

• number _propagat or s (integer) 

• rank (integer) 

• num(q,qt2) (complex function) 

• den0,denl ) den2,den3 ) den4,den5 (derived types: see below) 
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The first variable refers to the number of propagators in the (sub)-amplitude to be com- 
puted. The second variable is the maximum rank of N(q) (not greater than 
number_propagators, condition that is anyway always fulfilled in renormalizable gauges). 
num(q J qt2) is the numerator function N(q), that, when pieces of amplitude containing a 
different number of loop propagators are put together, also can depend on q 2 , that is the 
second entry of the function num(q,qt2). 

The last line of the above list refers to a derived type defined as follows 

module def _propagator 
implicit none 
type propagator 
integer : : i 
real (kind (l.dO)) : : m2 
real (kind(l . dO) ) , dimension(0 : 3) :: p 
end type propagator 
end module def _propagator 

Therefore, denj contains the information sufficient to denote the j th loop propagator, 
namely squared mass and 4-momentum. These loop propagators are internally classified 
according to a binary notation denj — > 2 J (following the user defined input ordering). 
The integer variable i of the previous derived type, is internally set to i = 2 s for each 
propagator. In the present version of CutTools, the maximum allowed number of loop 
propagators is six. When less propagators are needed, they should be loaded starting 
from the lowest value of j . 

At the end of the fitting procedure, the final results, namely the coefficients of the 
scalar loop functions and the rational part are loaded in the variables 

dcoeff (0, j) , ccoef f(0, j) , 

bcoeff (0, j) , bcoeff (3, j) , bcoeff (6, j) and rati. (16) 

The second index labels the relevant scalar loop functions, according to the above binary 
notation. For example the coefficient of the 3-point function 

^TO (17) 

is ccoeff(0, 2° + 2 2 + 2 4 ) = ccoeff(0, 21) and that one of 



/ 



d ^D^m < 18 > 

is dcoef f (0, 30). Furthermore, bcoef f (0, j), bcoeff (3, j), and bcoeff (6, j) are the coef- 
ficients of the scalar integrals in Eq. (14) with £ = 0, 1,2, respectively. When £ ^ 0, also 
the knowledge of the vector v is needed 4 . This information is stored in the array 

vvec(0:3,j), (19) 



4 The massless vector v is determined by CutTools in an event by event basis, to maximize the 
numerical stability. 
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where the second index j follows the same binary notation used for the loop propagators. 

Finally, when the multi-precision version of the code is activated, the relevant output 
information is stored in the variables: 

mp_dcoef f (0, j) , mp_ccoef f (0, j) , 

mp_bcoef f (0, j) , mp_bcoef f (3, j) , mp_bcoef f (6, j) , 

mp_vvec(0 : 3, j) and mp_ratl. (20) 

4 Program structure 

The directory structure looks as follows: 

avh_olo_s4.f dynamics. f 90 MPREC README 
cuttools . f 90 kinematics .f 90 process. f 90 tensors. f 90 
DOC Makefile rambo.f type. f 90 

In the following, we briefly discuss the content of each file or directory in the previous 
list. 

4.1 avh_olo_s4.f 

This set of routines, provided by Andre van Hameren, evaluates the scalar one-loop func- 
tions. In the current version the fully massless scalar one-loop functions are included [22]. 
However, when needed, since CutTools is not limited to massless processes, a more general 
repository of one-loop master integrals can be used [23]. 

4.2 cuttools.f90 

It is the main program. The distributed version implements, as a simple example, the 
reduction of a five-point function with a "toy" numerator 5 . A test run output is given in 
appendix B. 

The user should first initialize a few variables, such as the numbers of digits used by 
the multi-precision routines (idig), filling the internal tables of combinatorial factors by 
the calling the subroutine load_combinatorics, setting the the number of propagators 
for the case at hand (number_propagators), the maximum rank of N(q) (rank) and the 
limit of precision below which the multi-precision routines activate (limit). 

Then, for each generated Phase-Space point (the maximum number of points nitermax 
should be provided at running time), the user should define the derived types denj (j 
= 0, • • • , number_propagators-l) referring to the loop propagators, and load them by 

5 The routines for the evaluation of the complete one-loop QCD virtual corrections to the process 
qq — > ZZZ will also be available from our webpage. 
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calling the subroutine load_denominators (denO ,•••), with a number of arguments equal 
to the number of propagators. 

Finally, the needed coefficients of the one-loop scalar functions and the rational part 
Ri (rati) 6 are obtained by calling the subroutine get_coef f icients. At this point, 
if the precision test described in subsection 2.3 gives a result less then limit, the pro- 
gram multiplies all coefficients by the proper loop functions (this is achieved by calling 
dp_result (dbl_prec , cutpart) ), adds the rational parts and stores the event. Other- 
wise the entire procedure is repeated by using multi-precision. If the test fails even using 
multi-precision (that may happen if idig is too small), the event is discarded. 

At the end, the code, prints out the result of the Monte Carlo Phase-Space integration 
in the form of real and imaginary parts of the finite term (sigma(O)) and of the coefficients 
of the 1/e (sigma(l)) and 1/e 2 (sigma(2)) poles. A statistics is also provided of the 
percentage of points computed with multi-precision or discarded. 

4.3 DOC 

It is a directory containing this paper and any other updated documentation. 

4.4 dynamics. f90 

It is the part of the code where the user has to insert the numerator function N(q), namely 
the complex function num(q,qt2) . 

4.5 kinematics. f90 

It is the core of CutTools. It contains all routines needed to perform the fits. All the out- 
put variables listed in Eq. (16), Eq. (19) and Eq. (20) are located in module coefficients. 

4.6 Makefile 

It is the Makefile of CutTools. The user should specify, among other things, the 
F0RTRAN90 compiler and the compilation flags he/she is using. Notice that the multi- 
precision library in MPREC should be compiled first (see next subsection). 

4.7 MPREC 

It is a directory containing the multi-precision package of [21]. More precisely, before 
compiling CutTools, the user should go to /MPREC/mpf un90/f 90 and give the command 
make to compile the multi-precision library. 

6 We recall that the rational term i?2 (rat2) should be computed separately. 
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4.8 process. f90 

All routines needed to compute N(q) should be put in this file. 

4.9 rambo.f 

It contains the random number generator and the routines for Phase-Space generation, 
histogramming and bookkeeping of the events. 

4.10 README 

It is a .txt file with information on the current version of CutTools. 

4.11 tensors. f90 

It contains the routines needed to perform scalar products of 4-vectors. 

4.12 type.f90 

It contains the F0RTRAN90 derived types used by CutTools. 

5 Conclusion 

We have presented CutTools, a program implementing the OPP reduction method [15] 
to extract the coefficients of the one-loop scalar integrals from a user defined numerator 
function (namely (sub)-amplitude or Feynman Diagram), as well as the rational terms 
of type Ri [20]. The remaining part of the rational terms, R 2 , should be supplied by 
the user and can be computed with extra Feynman rules, as described in [20]. The 
possible occurring numerical instabilities are treated with the help of arbitrary precision 
routines [21]. The OPP algorithm allowed us to implement a trivial check in order to 
activate the time consuming arbitrary precision routines only when necessary. 
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Appendices 

A Going from double to multi-precision 

All routines in CutTools have been written both in a normal form (namely in double- 
precision) and in a multi-precision form. Once a routine is written in normal form, the 
multi-precision version of it can be easily obtained through the following changes in the 
declarations statements [21]: 

real(kind(l.dO)) — > type(mp_real) 

complex(kind(l.dO)) — > type(mp_ complex) . (21) 

The same strategy should be applied by the user to provide the multi-precision version of 
the routines to compute N(q). Finally, an interface statement can be used to call both 
versions with the same name. 

B Test run output 

With nitermax= 1, the final output of the program reads as follows: 
Result of the integration: 

real_sigma(0)= -67151075.5213172 +- 0.00000000000000 

imag_sigma(0)= -28426491.5346667 +- 0.00000000000000 

real_sigma(l)= 3822974.11389803 +- 0.00000000000000 

imag_sigma(l)= 3694493.63813566 +- 0.00000000000000 

real_sigma(2)= 1 . 925254707235981E-029 +- 0.00000000000000 
imag_sigma(2)= -1 . 427701757691647E-028 +- 0.00000000000000 

Statistics on the mp routines: 

percentage of mp points= 0.00000000000000 

percentage of discarded points= 0.00000000000000 

digits used in mp routines (if called) = 57 
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